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ABSTRACT 

A large proportion of observed planetary systems, which contain several planets in a compact orbital configuration, often harbor at 
least one close-in object. In this case, these systems are most likely tidally evolving. We investigate how the effects of planet-on-planet 
interactions influence the tidal evolution of planets. 

To achieve this, we introduced a new open-source addition to the Mercury N-body code, Mercury-T , which takes tides, general 
relativity (GR), and the effect of rotation-induced flattening into account to simulate the dynamical and tidal evolution of multiplanet 
systems. This code uses a standard equilibrium tidal model, the constant time lag model. Additionally, the evolution of the radius of 
several host bodies has been implemented (e.g., brown dwarfs, M dwarfs of mass 0.1 M©, Sun-like stars, and Jupiter). We validate 
the new code by comparing its output for one-planet systems to the secular equations results. We find that this code respects the 
conservation of total angular momentum. 

We then applied this new tool to the planetary system Kepler-62. As a result, we find that, in some cases, tides influence the stability of 
the system. We also show that, while the four inner planets of the systems are likely to have slow rotation rates and small obliquities, 
the fifth planet could have a fast rotation rate and a high obliquity. This means that the two habitable zone planets of this system, 
Kepler-62e and Kepler-62f, are likely to have very different climate features and of course, this influences their potential for hosting 
surface liquid water. 

Key words. Planets and satellites: dynamical evolution and stability - Planet-star interactions - Planets and satellites: terrestrial 
planets - System: Kepler-62 


1. Introduction 

More than 1400 exoplanets have now been detected and about 
20 % of them are part of multiplanet systems (http:// 
exoplanets.org/). Many of these systems are compact and 
host close-in planets where tides have an influence. In particu¬ 
lar, tides can have an effect on the eccentricities of planets, and 
also on their rotation periods and their obliquities, which are im¬ 
portant parameters for any climate study. Moreover, tides can 
influence the stability of multiplanet systems, as a result of their 
effect on both the planet’s eccentricities and precession rates. 

We present a new code, Mercury-T 1 , which is based on the 
N-body code Mercury (Chambers 1999). This allows us to cal¬ 
culate the evolution of semi-major axis, eccentricity, inclination, 
rotation period, and obliquity of planets, as well as the rotation 
period evolution of the host body. This code is flexible, in that it 
allows us to compute the tidal evolution of systems orbiting any 
non-evolving object (provided we know its mass, radius, dissi¬ 
pation factor, and rotation period), as well as evolving brown 


1 The link to this code and the manual can be found here: http: // 
www. emelinebolmont. com/. 


dwarfs (BDs), an evolving M dwarf of 0.1 M 0 , an evolving Sun¬ 
like star, and an evolving Jupiter. 

The dynamics of multiplanet systems with tidal dissipation 
have been the subject of study (evolution of the orbit in Wu & 
Goldreich 2002; Mardling 2007; Batygin et al. 2009; Mardling 
2010; Laskar et al. 2012 and also of the spin in Wu & Mur¬ 
ray 2003; Fabrycky & Tremaine 2007; Naoz et al. 2011; Correia 
et al. 2012), but most of these studies use averaging, do not study 
the influence of an evolving host body radius, and often consider 
only coplanar systems. In this paper, we introduce a tool that al¬ 
lows for more complete studies. Indeed, the tidal equations used 
in this code are not averaged equations, which makes it possible 
to study phenomena such as resonance crossing or capture. Con¬ 
trary to other codes using semi-averaged (Mardling & Lin 2002) 
or non-averaged equations (Touma & Wisdom 1998; Mardling 
& Lin 2002; Laskar et al. 2004; Fienga et al. 2008; Beauge & 
Nesvomy 2012; Correia & Robutel 2013; Makarov & Berghea 
2014; and Plavchan et al. 2015), our code is freely accessible and 
open source. 

After describing the tidal model used here, we seek to vali¬ 
date the code by comparing one planet’s evolutions around BDs 
with evolutions computed with a secular code (as in Bolmont 
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et al. 2011, 2012). We use systems around BDs to test systems 
where tides are very strong and lead to important orbital changes. 
We then offer a glimpse of possible research that could use our 
code and illustrate this with the example of the dynamical evo¬ 
lution of the Kepler-62 system (Borucki et al. 2013). 

2. Model description 

The major difference between Mercury and Mercury-T is the ad¬ 
dition of tidal forces and torques. However, we also add the ef¬ 
fect of general relativity and rotation-induced deformation. In 
the following sections we explain how these effects were incor¬ 
porated in the code. We also give the planets and star/BD/Jupiter 
parameters which are implemented in the code. 



2.1. Tidal model 

To compute the tidal interactions, we used the tidal force as ex¬ 
pressed in Mignard (1979), Hut (1981), Eggleton et al. (1998), 
and Leconte et al. (2010) for the constant time lag model. This 
model is based on the assumption that the bodies under review 
are made of a weakly viscous fluid (Alexander 1973). 

We added this force in the N-body code Mercury (Chambers 
1999). We also consider the tidal forces between the star and 
the planets but ignore the tidal interaction between planets. In 
addition, we consider here a population of N planets orbiting a 
star. 

As in Hut (1981), in order to obtain the expression of the 
force, we stop the development at the quadrupole order. At this 
order, we can use the point mass description of the tidal bulges. 
Star and planets are deformed. Owing to the presence of planet 
j, the star of mass M* is deformed and can be decomposed in 
a central mass M* - 2//*, and 2 bulges of mass p+. As in Hut 
(1981), each bulge is located at a radius 7?* from the center of 
the star and they are diametrically opposed. Figure 1 shows the 
geometrical context of the problem. The central mass of the star 
is labeled S, and the bulges S’ and S”. The mass of a bulge de¬ 
pends on the time lag and is given by 

F* = ( r j( f “ T *)) 3 . (!) 

where rj is the distance between the star and planet j at time 
t - r*, R+ is the radius of the star, k 2 ; * its potential Love number 
of degree 2, and r* its constant time lag. 

Because of the presence of the star, the planet j is deformed 
and can be decomposed in a central mass M Pj - 2fi V] , and two 
bulges of mass /i Pj . The central mass of the planet j is labeled by 
Pj, and the bulges by Pj ’ and Pj”. The bulge’s mass is given by 

^pj “ 2^ 2 ’ PjJ ^^^Pj ( r J^ ~~ Tp j)) ’ (2) 

where R Pj is the radius of planet j, & 2 ,pj its potential Love num¬ 
ber of degree 2, and r Pj its time lag. To the lowest order in r Pj , 
Equation 2 becomes: 


Fig. 1. Two-dimensional diagram representing the two deformed bod¬ 
ies. The star is divided into three masses: a central mass of M* - 2p+ at 
S, and two bulges of mass p+ at S’ and S”. The planet j is divided into 
three masses: a central mass of M p . - 2p v . at Pj, and two bulges of mass 
p V] at Pj' and Pj". £1* is the star rotation vector (its norm is Q*, the star 
rotation frequency), H p is planet j rotation vector (its norm is O p , the 
planet rotation frequency), and 6 } is a vector collinear with the orbital 
angular momentum of planet j (its norm is equal to the derivative of the 
true anomaly). e rj is the radial vector. 


where PjS" is the vector PjS", defined in Figure 1. 

Let us define 7q r (for tides radial), P to> *, and Z\o,pj (f° r tides 
ortho-radial) as 

F tt = {M 2 p .h,*Rl + Mlk 2 , p Rl) 

J 

- 90 V -F + MlR 5 p k 2 , Pi T Pi ), 

J 

MlRl 

^Vpj = 3£r k 2 , pjTpj, 

J 

MlRl 

P t0 ,* = 3g-?L—*2,*r*. (5) 

r! 

j 

Here, F a has the dimension of a force (M.L.T -2 ), while P t o,pj and 
P t0 ,* have a dimension of a momentum (M.L.T -1 ). 

Consequently, the total resulting force as a result of the tides 
acting on planet j is 


Fa + [P to,* 

+ Ft o,pj (A, 
+ P to,* (f^* 


0j) x e rj 
- 0j) x e rj , 


(6) 






V 3 


1 + 3— Tn 


(3) 


Up to the third order in R p Jr s and R./i], the forces exerted 
by the primary on the secondary are the following gravitational 
forces: f s ^ Pj , f's^ps fs->p", fs'^p,. and f S "^p,, where the latter 
expression is given by 


(iUp, 2 A*Pj) p 

l|PjS"|| 3 j 


(4) 


where O* is the star rotation vector, Q Pj is planet j rotation vec¬ 
tor, and Vj = rj is the velocity of planet j. The unit vector, e Fj , is 
defined as SPj/SPj, while is a vector collinear with the orbital 
angular momentum of planet j (defined hereafter as LhorbX the 
norm of which is equal to the derivative of the true anomaly. The 
term 0j x e rj can be re-written as follows: 

6\ x e,. = f (rj x vj) x e rj . (7) 
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What modifies the spin of the star is the following torque 
contribution: -r j X (f,^.s- + fpj-S"). To calculate the torque on 
the star, N p .^*, we consider the planet as a point mass, meaning 
that we ignore the gravitational interaction between the bulges of 
the planet and the bulges of the star. Following a similar hypoth¬ 
esis, we find that the torque contribution that modifies the spin 
of planet j, N*_, p ., is rj x so the torque exerted 

by planet j on the star is given by 

N£_* = N* = P to ,* (rj fi* - (rj.il*) e rj - e rj x Vj), (8) 

and the torque exerted by the star on the planet j is 

= Np. = Pto.pj ( r j ^pj “ ( r j *^Pj) e rj “ e rj XVj). (9) 

With this description of the phenomenon, we consider that 
each planet creates an independent tidal bulge on the star and 
that the bulge created by planet j does not affect planet iAj. 


2.2. General relativity 


We added the force due to general relativity as given in Kid¬ 
der (1995); Mardling & Lin (2002) to Mercury. This force cor¬ 
responds to the orbital acceleration due to the Post-Newtonian 
potential and its expression is 

f£ r = M Pj (F GRl .e rj + F GRo vj). (10) 

Here, Fgr, and Fqro are given by 


^GRr = -- 


£(M* + M Pi ) 


X ^1+377)^-2(2 + 77)- 

n _ 0/0 G(M+ + M Pj ) 

^gro - 2(2 - rj) -—--rj, 

r L c l 

J 


£(M* + M Pi ) 3 




(ii) 


where vj is the norm of the velocity Vj of the planet, and c is the 
speed of light, and where rj is 


m*m Pj 

(M* + M Pj ) 2 - 


( 12 ) 


2.3. Rotational deformation 


The equilibrium figure of a viscous body in rotation is a triaxial 
ellipsoid symmetric with respect to the rotation axis (Murray & 
Dermott 1999). The rotational deformation is quantified by the 
parameter J 2 , defined for planet j as follows : 


A Pl - % P| ■ 

and for the star as 


(13) 


J 2 ,* - kif,* 


njRi 

3 £M*' 


(14) 


Here, £ 2 /^ is the fluid Love number of planet j and & 2 /> that of 
the star. We define the fluid Love number as the potential Love 
number for a perfectly fluid planet (see, for example, Figure 2 of 
Correia & Rodriguez 2013, for the Earth’s potential Love num¬ 
ber and fluid Love number). Our code allows the user to choose 


different values for the fluid Love number & 2 /,p and the potential 
Love number & 2 ,p- 

The total resulting force due to the rotational deformation of 
star and planet j on planet j is (Murray & Dermott 1999; Correia 
etal. 2011) 


j 


(c* + c Pj ) 

rj.il* 


15 


C* 


M*f . „ Mr) 




ft* 


+ Crv 


o 2 


rj.n D 

n * + c p j -^r LQ P 


where C* and C Pi are defined as follows: 


1 


C+ — ^ /2,pj ^ Pj 


c Pj = ] -gM p M+J 2 _.Ri. 




N 


The torque exerted by planet j on the star is given by 

6_ r r j ft* 

p i-* "* r 5 C * q2 


n; 


(rj x il*), 


(15) 


(16) 


(17) 


and the torque exerted by the star on the planet is 


N r 


- n r - --C 

- iA Pj - r 5 C Pj 


r l -ftp, 
Or 

^Pj 


(rjxftpj. 


(18) 


These torques are responsible for the precession of the orbit 
normal. This precession has an influence on the mean eccentric¬ 
ity of planets and also on their mean obliquity. 


2.4. Summary of all the effects 
2.4.1. Corrective acceleration 

To compute the evolution of the planetary orbits, we need to take 
all the resulting effects into account. The orbital part of the accel¬ 
eration is handled by Mercury , so below we provide the expres¬ 
sion of the corrective acceleration of planet j in the astrocentric 
coordinates: 


M Pj + M* 

a D . =- 

Pj M Pj M* 


K + f b r+ r“) + ^ 2 (fJ+fS r +f;) 

* i*J 


_L ( f t + f gr + f r 1 + _L y ( f t + ¥ gr f r\ 

M v Pj Pj Pj/ M V Pi Pi Pi/ ’ 

iKZ Pj 1V1 * 7=1 


(19) 


where F^, F p ’ R and F R are defined in the previous sections. 


2.4.2. Spin equations 


First, let us consider a system with one planet. We hypothesize 
that we can decouple the torque equation given by the conserva¬ 
tion of total angular momentum L. This equation is the follow¬ 
ing: 


d t 
d 
d t 


L = 0, 


+ L 0 rb) - 0, 


( 20 ) 
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where I * is the principal moment of inertia of the star, 7 p that of 
the planet and L or b is the orbital angular momentum: 


Lorb = r* x MpV* + < x M*v;, 


( 21 ) 


where r* and v* are the position and velocity of the planet in the 
reference frame of the center of mass of the system. The position 
and velocity of the star in the reference frame of the center of 
mass of the system are respectively r* and v*. 

In astrocentric coordinates, the position and velocity of the 
planet are r and v and Equation 21 becomes: 


L 


orb — 


M*M p 

--rxv, 

M* + M p 


( 22 ) 


so 


d 

dr 



M*M P d I 

-r rxv 

A7, + Alp dr v . 


M* 

M* + M p 


(N*-> p + N p _>*), 


(23) 


where N+^ p is the total torque exerted on the planet, and N p ^+ 
the total torque exerted on the star. 

Assuming that the spins of the star and the planet evolve 
solely as a result of tidal and rotational flattening torques, this 
means that we can decouple Equation 20 to obtain the following 
spin equations: 


&('*«*) 


d_ 

d t 



=-iot( n ; +n ?). 


(24) 


If there is more than one planet in the system, the orbital an¬ 
gular momentum involves planet j- planet iAj cross terms (e.g., 
in Fj x Vj, Fi x Vj) that we ignore for the spin calculation. The 
equations governing the evolution of the spin of the star and the 
spin of planet j are, therefore, 

£('*«*) =-2s^(nI + n5) 

J= 1 P) (25) 

>('«,“») =-»^K + N E)' 


where N*, N*, N p ., and are defined in the previous sections. 
The spin of planet j is the norm of the vector Q Pj and the spin of 
the star is the norm of the vector O*. 

Our Mercury-T code provides the x, y, and z components of 
the spin of the planets, of their orbital angular momentum, and 
of the spin of the star. Then the obliquity £ Pj and the inclination 
/j of planet j are obtained by calculating 


cos e n; 


Lorbj ' ^pj 


-'orb; 


x £2 n 


COS /j 


Lorbj ' A* 

IlLorlJ X ||n*|f 


(26) 


where is a vector normal to the orbit of planet j. 
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2.5. Integration of the spin 

In this study, we use the Mercury code’s hybrid routine that relies 
on a Hamiltonian description of the problem. The Hamiltonian is 
divided into three parts that are integrated consecutively (Cham¬ 
bers 1999). Mercury allows the user to add other forces on the 
planets via a routine. In the Hamiltonian description, these ex¬ 
tra forces are treated as a perturbation to the Keplerian potential. 
The user should keep in mind that this violates the symplectic 
properties of the integrator. 

In this user routine, we add the tidal forces, rotation-induced 
flattening forces and GR forces. On the one hand, we computed 
the resulting acceleration on the planets, and this acceleration 
is used by the bulk of Mercury to compute the evolution of the 
system. On the other hand, the spin of planets and the star is 
integrated with a 5th order Runge-Kutta within the routine. The 
Runge-Kutta integration is performed twice in a Mercury time 
step. To compute the spins at a time t , the positions and velocities 
of the planets are interpolated between t-dt and t at the time 
intervals required for the Runge-Kutta routine. The integration 
scheme is illustrated in Figure 2. 


t-dt 


t 


Mercury 

x^(t-dt) 

vj(t-dt) 


3p,tot(t-dt) 


Mercury 

Xp(t) 

Vp(t) 




routine 

k. — 

Qp(t-dt) 

Np 

n P (t) 

n*(t-dt) j\T* 

sS*(t) 


Fig. 2. Integration scheme of Mercury-T. 


If the host body evolves, the moment of inertia of the star, 
given by: I * = M+(rg+R+) 2 , where rg+ is the radius of gyration 
(Hut 1981), varies with time. 

The equation of each component of the spin O* is given in 
the following equation, written here for the z component 


/*(0^*, z (0 = hit - d t)£l+, z (t - dt) 



M* + M Pj 


{Kz + Kz)dt- 


(27) 


The error introduced by the spin integration depends on 
the third power of the time step (Chambers 1999). The part of 
the integration that causes more errors is the integration of the 
rotation-induced flattening effect. The integration of the tidal 
torque does not require such a precise integration owing to long 
timescales of evolution with respect to the Mercury time step 
(which is usually taken as slightly smaller than one tenth of 
the inner planet’s orbital period). However, the rotation-induced 
flattening causes changes in the spin of the planets on a much 
shorter timescale. 

For very close-in planets, the integration of the rotation- 
induced flattening effect can lead to a purely numerical decrease 
in the rotation period. It is therefore important to evaluate the 
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error made during the integration of the rotation-induced flatten¬ 
ing torque. We suggest using a time step shorter than the orbital 
period of the inner planet divided by 20. The time step should 
therefore be chosen according to the precision required in the 
rotation period of the inner planet (see section 3.5 for details). 

2.6. Input parameters 

The code requires the necessary planetary parameters to work 
and these are used in all the equations in sections 2.1, 2.2, 
and 2.3. An N-body integrator requires parameters such as: the 
masses of the planets M Pj , the mass of the host body M+, the 
semi-major axis (SMA) of the planets, their eccentricities (ecc), 
their inclination (inc), and their orbital angles (argument of peri- 
center, longitude of the ascending node, and mean anomaly). In 
the following tests, the orbital angles are set to 0°. 

To calculate the evolution that results from rotational flatten¬ 
ing, the fluid Love number of the star k 2 /> and of the planets 
k 2 , Pj is required. To calculate the tidal evolution, the following 
are also required: the radius of the star R* and the radius of the 
planets R Pj ; the potential Love number of degree two of the star 

and of the planets k 2 , Pj ; the time lag of the star r* and of the 
planets r Pj . 

All these parameters are, of course, changeable in our code, 
but we implemented some useful values and relations for ease of 
use. These implemented values are given in Tables 1 and 2. 

2.6.1. Planet model 

For Earth-like planets and super-Earths, we assume that the 
product of the potential Love number of degree two with the time 
lag of the planet is equal to that of Earth: k 2 , p Ar p = k 2 ,©Ar®. We 
use here the value of k 2 ,©AT© = 213 s given by Neron de Surgy 
& Laskar (1997). We assume here that the fluid Love number 
and the potential Love number of degree two are equal. 

Given the mass of the planet, we offer the user two possibili¬ 
ties to choose the radius: either it gives a value itself or it assumes 
a composition and the code calculates the radius following Fort¬ 
ney et al. (2007). For example, a super-Earth of 10 M© would 
have a radius of 1.8 R©. 

For a Jupiter-like planet, we computed the time lag thj from 
the value of the dissipation parameter <Xk for hot Jupiters of 
Hansen (2010). The notation <Xk was introduced by Eggleton 
et al. (1998) and is linked to the quantity /^k^k hy 


^2,k^"k — 


3 fljVk 
2 ~Q~' 


where k represents either the star or a planet. 


(28) 


2.6.2. Host body evolution and dissipation 

It is possible to use stellar evolution tracks in Mercury-T to 
compute the evolution of planets around an evolving object. We 
implemented this for an evolving Jupiter (Leconte & Chabrier 
2013), and for evolving BDs of masses: 0.01, 0.012, 0.015, 0.02, 
0.03, 0.04, 0.05, 0.06, 0.07, 0.072, 0.075, and 0.08 M© (Leconte 
et al. 2011), for a M dwarf (dM) of mass 0.1 M© and a Sun-like 
star (Chabrier & Baraffe 1997; Baraffe et al. 1998). 

Table 2 shows for each type of host body which ones are the 
evolving quantities and which ones have implemented values. 
The evolving quantities are tabulated and Mercury-T interpolates 
the values during the integration to have the correct radius, Love 


number, moment of inertia in the acceleration formula and the 
spin equations. 

At this point, we hypothesize that, during their evolution the 
dissipation factor, cr k of the host body remains constant. We use 
the value of the dissipation of Jupiter given in Leconte et al. 
(2010) for the Jupiter host body. We use the dissipation factor 
of Bolmont et al. (2011) for BDs and for the M dwarf, and we 
use the value of stellar dissipation of Hansen (2010) for the dis¬ 
sipation of the Sun-like star. 

In Table 3, we indicate the value of the parameter erf for the 
different bodies, where erf is defined as 


ov = M.RiP n <T+, (29) 


and where = ^jR* IQM+ is the free-fall time at the surface 

of the star. As this definition of erf depends on the radius of the 
star, we compute its value for all the bodies of Table 2 for an age 
of 1 Myr and of 1 Gyr. 


3. Code verification 

To validate the tidal part of the code, we first simulated the tidal 
evolution of one Earth-mass planet orbiting a 0.08 M© BD with 
two different approaches. The first approach is to use a secular 
code that solves the averaged equations of the tidal evolution 
of one planet (equations in semi-major axis, eccentricity, etc, 
which is often used in tidal studies such as Hut 1981; Leconte 
et al. 2010; Bolmont et al. 2011). The second simulation was per¬ 
formed using the Mercury-T code we developed. We first com¬ 
pared the outcomes for different regimes, switching the planetary 
tide on or off (i.e., the tide raised by the BD in the planet) and 
the BD tide (i.e., the tide raised by the planet in the BD), and 
testing the effect of the evolving radius of the BD. The details 
of all simulations are listed in Table 4, where d t is the time step 
used for the simulation. 

To test the rotational flattening part of the code, we compared 
our results with the numerical code used in Correia & Robutel 
(2013), hereafter denoted by the CR13 code. This code was de¬ 
veloped independently of the present one and uses the ODEX 
integrator (e.g., Hairer et al. 1993), but it has not been made 
available for public use. In Correia & Robutel (2013), the CR13 
code was applied to a specific situation: the spin evolution of 
trojan bodies, but it is much more general than that. The CR13 
code has the ability to perform the same kind of simulations as 
those of the Mercury-T code. Therefore, the CR13 code is used 
for cross-checking some of our results. 

We considered a two-planet system orbiting a 0.08 M© BD 
and validated the effect of the rotational flattening of the star, 
and of the planet and compared them with our results for a full 
simulation (effects of tides and rotational flattening). The details 
of the simulations are listed in Table 5, where & 2 , Pl/2 is the Love 
number of degree 2 of planets 1 and 2, and is the Love num¬ 
ber of degree 2 of the host body. 

3.1. Non-evolving BD: the effect of planetary tide 

This case corresponds to Case (1) of Table 4, for which we 
switched off the effect of the BD tide. 

Figure 3 shows the evolution of the semi-major axis, the ec¬ 
centricity, and the averaged tidal heat flux (0 t ides) defined as 

(0tides) — <£tides>/47TR p , (30) 
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Table 1. Planetary parameters implemented in Mercury-T. 


Type of planet 

Mass 

Radius 

Love 

number 

Moment of inertia 

I/(MR 2 ) 

Time 1 
Value (s) 

ag 

Notation 

Earth-mass planet 
Jupiter 

1 M e 

lMj 

I*® 

iflj 

0.305 

0.380 

0.3308 

0.254 

698 

1.842 x 10“ 3 

T© 

THJ 


Table 2. Host body parameters implemented in Mercury-T. 


Type of 

Mass 

Radius 

Love 

Moment of inertia 

Dissipation factor 

host body 



number 

// (MR 2 ) 

Value (g l cm 2 s l ) 

Notation 

Jupiter 

lMj 

evolving 

evolving 

evolving 

7.024 x 10“ sy 

crj 

BD 

0.01-0.08 M e 

evolving 

0.379-0.307 

evolving 

2.006 x 10- 60 

O' BD 

dM 

0.1 M 0 

evolving 

0.307 

0.2 

2.006 x 10“ 60 

O' d M 

Sun 

1M g 

evolving 

0.03 

0.059 

4.992 x 10“ 66 

o~o 


Table 3. Values of the reduced dissipation factor cr* for the host bodies implemented in Mercury-T. 


Type of 

Mass 

Radius 


host body 


at t = 1 Myr 

U 
>, 
o 

II 

C3 

c\it= 1 Myr 

at t = 1 Gyr 

Jupiter 

lMj 

0.15 Rq 

0.10 Rq 

4.3 x 10 -5 

1.1 x 10- 3 

BD 

0.01 M q 

0.29 Rq 

0.10 Rq 

4.0 x 10“ 5 

9.8 x 10“ 7 


0.08 Mq 

0.85 Rq 

0.10 Rq 

4.9 x 10“ 3 

2.7 x lO -6 

dM 

0.1 Mq 

0.98 Rq 

0.12 Rq 

9.0 x 10“ 3 

5.8 x 10“ 6 

Sun 

1 Mq 

2.3 Rq 

0.91 Rq 

1.4 x 10 -6 

5.5 x 10“ 8 


a) b) 



Time (years) Time (years) 

Fig. 3. Case (1): Tidal evolution of a planet of mass 1 M© orbiting a 0.08 M© BD, calculated with the secular code (blue dashed line) and 
the Mercury-T code (solid red line). Graph a) from top to bottom: evolution of the semi-major axis of the planet, evolution of its eccentricity, 
and evolution of the tidal heat flux. Graph b) from top to bottom: evolution of the obliquity of the planet, evolution of its rotation period (the 
pseudo-synchronization period represented in long red dashes), and conservation of total angular momentum. 


where (E t id es ) is the averaged gravitational energy lost by the 
system by dissipation. Here, 


(^tides) —2 


1 QM V M+ 


4 a 


D p 

Nal(e) - 2Na2(e) cos 6 P — 
n 


' 1 + cos 2 e p \ /a,' 2 

+ , ^7Tl Q(e) f 
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( 31 ) 


where T v is the dissipation timescale. Nal(e), Na2(e ), and £l(e) 
are eccentricity-dependent factors defined in Bolmont et al. 
(2013). The tidal heat flux depends on the eccentricity and on 
the obliquity of the planet. If the planet has no obliquity and no 
eccentricity and if its rotation is synchronized, the tidal heat flux 
is zero. 
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Table 4. Test simulations for tides: one BD and one planet. 




Effects 

Parameters and initial conditions 


Tides 

BD 

M* 

R+ 

p*. 0 

M p 

SMA 

ecc 

inc 

Pp 

e p 

t p 

<r* 

d t 


BD 

PL 

evol. 

(M 0 ) 

(Re) 

(day) 

(M e ) 

(AU) 


(deg) 

(hr) 

(deg) 



(day) 

1 

X 

/ 

X 

0.08 

- 

- 

1 

0.014 

0.1 

0 

24 

11.5 

T© 

- 

0.08 

V 

X 

/ 

X 

0.08 

~ 

- 

318 

0.014 

0.01 

0 

24 

11.5 

100 r H j 

- 

0.08 

2 

/ 

X 

X 

0.08 

0.85 

2.9 

1 

0.018 

0 

11.5 

- 

- 

- 

1000 (T BD 

0.08 

3 

/ 

/ 

X 

0.08 

0.85 

2.9 

1 

0.018 

0.1 

5 

24 

11.5 

T® 

& BD 

0.08 

3’ 

/ 

/ 

X 

0.08 

0.85 

2.9 

318 

0.009 

0.1 

5 

240 

40 

100 r H j 

0.01 cr BD 

0.005 

4 

/ 

/ 

/ 

0.08 

evol 

2.9 

1 

0.018 

0.1 

0 

24 

11.5 

T© 

c BD 

0.05 


Table 5. Test simulation for rotational flattening: one BD with 2 planets 



Effects 

Parameters and initial conditions for planets 1 and 2 


Rot. flat. 

M* 

^Pl/2 

sma 1/2 

ecci/2 

inci/2 

^Pl/2 

6 Pl/2 

hpin 


d t 


BD 

Planets 

(M g ) 

(M e ) 

(AU) 


(deg) 

(hr) 

(deg) 



(day) 

5 

/ 

X 

0.08 

1/1 

0.018/0.025 

0.01/0.01 

o/i 

24/24 

11.459/11.459 

- 

0.307 

0.05 

6 

X 

/ 

0.08 

1/1 

0.018/0.025 

0.01/0.01 

0/1 

24/24 

11.459/11.459 

0.305/0 

- 

0.08 

6’ 

X 

/ 

0.08 

1/1 

0.018/0.025 

0.01/0.01 

o/i 

24/24 

11.459/11.459 

0.305/0 

- 

0.05 

6” 

X 

/ 

0.08 

1/1 

0.018/0.025 

0.01/0.01 

o/i 

24/24 

11.459/11.459 

0.305/0 

- 

0.01 

6’” 

X 

/ 

0.08 

1/1 

0.018/0.025 

0.01/0.01 

0/1 

24/24 

11.459/11.459 

0.305/0 

- 

0.001 

7 

/ 

/ 

0.08 

1/1 

0.018/0.025 

0.01/0.01 

o/i 

24/24 

11.459/23 

0.305/0.305 

0.307 

0.08 


We also show the evolution of the instantaneous tidal heat 
flux, computed from the instantaneous energy loss given by 

£tides(0 = -£orb(0 = - ( F r ' v j + 4j n Pj ' W ) ’ (32) 

where Q Pj is the derivative of the spin of planet j, given by Equa¬ 
tion 25. Contrary to (Elides) which depends on averaged com¬ 
puted values, such as the semi-major axis and the eccentricity, 
Elides(0 depends on the instantaneous position, velocity, and spin 
of planet j. 

The eccentricity of the planet decreases to values below 10~ 4 
in 10 7 yr. The decrease of the eccentricity of the planet is accom¬ 
panied by a decrease of the semi-major axis. The evolution of 
these two quantities shows a good agreement between the sec¬ 
ular code and the Mercury-T code. The obliquity of the planet 
decreases from its initial value of 11.5° to less than 10 -4 de¬ 
grees in less than 500 yr. During the same time, the rotation 
period evolves from its initial value of 24 hr to the pseudo¬ 
synchronization period, which in this study is ~ 48.5 hr. The 
evolution of obliquity and rotation period show a good agree¬ 
ment between the secular code and Mercury-T. 

After 2 x 10 7 yr of evolution, the eccentricity obtained with 
Mercury-T is equal to a few 10 -7 . This residual value of the ec¬ 
centricity comes from the way the Mercury code calculates the 
orbital elements from the positions and velocities of the planets. 
Indeed, it assumes a Keplerian potential. However, in this situa¬ 
tion, where the tidal forces are taken into account, this is not the 
case. 

Nevertheless, we can assume that an eccentricity of 1CT 7 be 
considered as null. Furthermore, this code is designed to study 
multiplanet systems. In the examples we give later, the eccen¬ 
tricity due to planet-planet interactions is typically greater than 
1(T 7 . 

This residual eccentricity is responsible for a non-zero- 
averaged tidal heat flux < 10 -2 W/m 2 . However, the instanta¬ 
neous tidal heat flux reaches values as low as a few 10 -9 W/m 2 . 


This low value illustrates the fact the real eccentricity of the 
planet must be much lower than that which Mercury-T is cal¬ 
culating. 

We also verify that each component of the total angular mo¬ 
mentum is a conserved quantity during the evolution of the sys¬ 
tem (Equation 20). We define the quantity a t where i is x, y, or 
z, and a as 


<*i = 


U(t)-U(0) 

L(0) 

L(0-L(Q) 
L(0) ’ 


(33) 


where L* is the i component of the total angular momentum vec¬ 
tor and L is the norm of the vector L of Equation 20. In this 
example, we only consider the effect of the planetary tide, which 
is equivalent to considering the BD as a point mass, so that the 
total angular momentum, in this case, is only the sum of the or¬ 
bital angular momentum and the rotational angular momentum 
of the planet. 

The bottom right panel of Figure 3 shows the conservation 
of the total angular momentum as a function of time. For the 
Mercury-T simulation, each component of the total angular mo¬ 
mentum oci is conserved and a reaches 10 -6 after 100 Myr of 
evolution. For the secular code, the total angular momentum is 
a little less well conserved and reaches a few 10 -5 at the end of 
the simulation. In this example, the orbital angular momentum 
is 10 5 higher than the angular momentum of the planet, so the 
question remains as to whether the spin of the planet has been 
correctly computed. 

To test this, we performed another simulation, Case (1’), with 
a Jupiter-mass planet to reduce the difference between the orbital 
angular momentum and the planet’s angular momentum. In this 
case, the orbital angular momentum is about 10 3 higher than the 
angular momentum of the planet. We find that for this simulation 
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Fig. 4. Case (3): Tidal evolution of a planet of mass 1 M© orbiting a 0.08 M© BD, calculated with the secular code (blue dashed line) and the 
Mercury-T code (solid red line). Graph a) from top to bottom: evolution of semi-major axis (in red) and of the corotation distance (in black), 
evolution of eccentricity, and evolution of tidal heat flux. Graph b) from top to bottom: evolution of the obliquity of the planet, evolution of its 
inclination and evolution of its rotation period (in red) and the BD rotation period (in black). The pseudo-synchronization period is represented in 
long red dashes. 


the total angular momentum is conserved, with a asymptopically 
reaching only 6 x 10 -6 after a 10 Myr evolution 2 . 

3.2. Non-evolving BD: the effect of BD tide 

This case corresponds to Case (2) of Table 4, where we switched 
off the effect of the planetary tide and considered a very dissipa¬ 
tive BD. Initially, the planet is outside the corotation radius, so it 
migrates outward. 

In agreement with the secular code, the BD tide causes the 
inclination of the planet to decrease from ~ 12° to ~ 4.5° in 
10 8 yr. As the planet migrates outward, the rotation period of the 
BD increases in agreement with the conservation of total angular 
momentum. Indeed, we find that every component of the total 
angular momentum is conserved. In this example, the angular 
momentum of the BD is of two orders of magnitude higher than 
the orbital angular momentum of the planet. Because a remains 
below 10 -4 after 100 Myr, we can conclude that the phenomenon 
is accurately reproduced here. 

3.3. Non-evolving BD: the effect of both tides 

This case corresponds to Case (3) of Table 4. In this example, 
the planet is initially outside the corotation radius. 

Figure 4 shows the evolution of this system. We find that 
the results of Mercury-T agree well with those using the secular 
code. The evolution of the different quantities are similar and 
the quantitative agreement is very good. As the planet migrates 
away, the evolution timescales become longer, which entails a 


2 On our server, the computation of this case, with a time step of 

0.08 day, required about five days to reach 10 Myr. This time would, of 

course, increase with more than one planet in the system and it would 

probably change on another computer. 


slower evolution, particularly in the late ages of the eccentricity 
and inclination. 

Moreover, the initial heat flux is very strong in compari¬ 
son to the tidal heat fluxes measured for solar system bodies: 
0.08 W/m 2 for Earth (Pollack et al. 1993) and between 2.4 and 
4.8 W/m 2 * * for Io (Spencer et al. 2000). Such a large heat flux 
is likely to have repercussions on the planet’s internal structure. 
The high fluxes of Figure 4 suggest that the surface and the in¬ 
terior of the planet would melt and that the vertical heat transfer 
could be very efficient, which does not agree with the dissipa¬ 
tion factor value used here. As in Bolmont et al. (2013), we do 
not include in this work any feedback of the dissipation on the 
internal structure. 

For this system, each components of the total angular mo¬ 
mentum ai is conserved and a reaches a few 10“ 7 at a time of 
10 Myr for the Mercury-T simulation. The total angular momen¬ 
tum here is, therefore, also conserved. 

We also test the strength of our code with a more extreme 
case, for example, Case (3 ? ). With an initial orbital distance of 
9 x 10 -3 AU, the planet is initially inside the corotation radius 
and thus migrates inward. Besides, we consider a BD of a low 
dissipation factor (0.01 xcr BD ) so that the planet’s BD-tide driven 
inward migration is not too quick. 

Figure 5 shows the evolution of this system. The planet 
plunges onto the BD in about 3000 yr. During the inward migra¬ 
tion, the eccentricity, obliquity, and inclination all decrease. In 
less than 2000 yr, the rotation period of the planet evolves from 
240 hr to the pseudo-synchronization period (of about 24 hr). 
The rotation period of the BD decreases just prior to the fall as a 
result of the angular momentum transfer from the planet’s orbit 
to the BD spin. 

Even for this extreme case, both codes lead to the same sim¬ 
ulated evolution. The collision time may be slightly different, 
but is of the same order of magnitude for both simulations. The 
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Fig. 5. Case (3’): Tidal evolution of a Jupiter-mass planet orbiting a 0.08 M 0 BD, calculated with the secular code (blue dashed line) and the 
Mercury-T code (solid red line). Graph a) from top to bottom: evolution of semi-major axis (in red) and of the corotation distance (in black), 
evolution of eccentricity and conservation of angular momentum. Graph b) from top to bottom: evolution of the obliquity of the planet, evolution 
of its inclination and evolution of its rotation period (in red) and the BD rotation period (in black). The pseudo-synchronization period is represented 
in long red dashes. 


bottom-left panel of Figure 5 shows the conservation of total an¬ 
gular momentum for this example. 

Because the planet ends up colliding with the BD, we do not 
expect the conservation of total angular momentum to be perfect. 
Indeed, Figure 5 shows that for both simulations, a increases 
with time and reaches about 10 -2 when the collision occurs (4 x 
10 -3 for Mercury-T). The planet is initially very close to the BD 
and gets closer in time, meaning the tidal effects become stronger 
and stronger. This example enables us to test the limits of our 
model. For close-in planets, one should always verify that a is 
conserved. 

In the end, the destiny of the planet is compatible with the 
theory. As its initial orbital distance is less than the corotation 
distance, the BD tide acts to push the planet inward. The qualita¬ 
tive evolution is not likely to change even if the code were to be 
improved, however the time of collision between the planet and 
the BD might change. 

3.4. Evolving BD: the effect of both tides 

In Case (4) of Table 4, we consider the evolution of the radius 
and radius of gyration of the BD. 

Our code allows us to choose the initial time of the simula¬ 
tions, i.e., the BD age from which we consider the tidal evolution 
of the planets. In this study, we assume that the initial time cor¬ 
responds to the time of the dispersal of the gas protoplanetary 
disk (as in Bolmont et al. 2011, where they discuss the influ¬ 
ence of this initial time). We then consider that the planets are 
fully formed by this time. The time indicated in the figures cor¬ 
responds to the time spent after this initial time. 

Figure 6 shows the evolution of this system. The evolution 
calculated with the secular code is in good agreement with the 
evolution calculated with Mercury-T . The competition between 


the outward migration caused by the BD tide and the inward 
migration caused by the planetary tide, is well reproduced. 

However, the comparison between the outcomes of the two 
codes shows a difference of 10 -5 AU in the calculated semi¬ 
major axis, when the migration direction changes. This differ¬ 
ence remains small and tends to decrease when the precision of 
the secular code is increased, which again demonstrates that the 
the Mercury-T code seems more precise than the secular code. 

In any case, the qualitative behavior is reproduced very well, 
even though small quantitative differences can be seen. The 
Mercury-T code reproduces the evolution of the spin of the BD 
well because of the contraction of its radius (middle panel of 
Graph b) in Figure 6). Besides, the total angular momentum is 
well conserved as can be seen in Figure 6. Indeed, each compo¬ 
nent of the total angular momentum as well as a remain below 
3 x 10“ 6 . 

These diverse tests show that the tidal integration part of the 
Mercury-T code shows a good agreement with the secular code 
in relation to the orbital evolution of the planet, as well as its ro¬ 
tation state evolution and the rotation evolution of the BD. The 
total angular momentum is always conserved, except when the 
planet collides with the BD. We therefore consider that this code 
is valid when studying the evolution of tidally evolving multi¬ 
planet systems. For any simulation, however, one should always 
make sure that the total angular momentum is conserved. 

3.5. Effect of the rotational induced flattening 

While an Euler integration of the spin may have been sufficient 
to correctly describe the tidal evolution of the spin of plan¬ 
ets, we need to implement a better integrator to accurately de¬ 
scribe the precession of the planet’s spin axis that results from 
its own flattening. This precession happens on a much shorter 
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Fig. 6. Case (4): Tidal evolution of an Earth-mass planet orbiting a 0.08 M 0 BD, calculated with the secular code (blue dashed line) and the 
Mercury-T code (red line). Graph a) from top to bottom: evolution of semi-major axis (in red) and of the corotation distance (in black), evolution 
of eccentricity, and evolution of tidal heat flux. Graph b) from top to bottom: evolution of the obliquity of the planet, evolution of its rotation period 
(in red), the BD rotation period (in black), and the pseudo-synchronization period (red dashed line), and evolution of a. 


timescale than the tidal evolution. For the example, in Case (5) 
the timescales are about a few 10 1 yr. 

Thus, to obtain an accurate integration with Mercury-T , we 
need to perform a 5th order Runge-Kutta integration twice in 
a Mercury time step. Dividing the time step into two inside one 
Mercury time step allows us to be more precise, without demand¬ 
ing too much time. 

We tested the integration of the rotation-induced flattening 
by comparing our code to the CR13 code. In doing so, we obtain 
similar results for all cases, with a few quantitative differences. 

For Case (6) of Table 5, only the rotational flattening of the 
inner planet is taken into account. The rotation period P v (i.e., 
the norm of the spin) is not influenced by the effect of the ro¬ 
tational flattening. However due to the integration scheme, we 
observe a small drift in the rotation period of the inner planet 
(Figure 7). This drift increases linearly with time and decreases 
when the time step is reduced. For a time step of 0.08 day, i.e., 
Case (6), the drift is of about 8 x 10 -5 hr after 100 000 yr of 
evolution, while, for a time step of 0.01 day, i.e., case (6”), it is 
less than 10 -6 hr. For time steps shorter than 0.01, the drift is 
essentially null for 100 000 yr of evolution. 

The shorter the time step, the smaller the differences. From a 
time step of 0.08 day to 0.01 day the improvement is visible, but 
we can see that there is almost no difference between the light 
and dark blue curves corresponding to time steps of 0.01 and 
0.001 day. Of course, the execution time is longer for shorter 
time step, so the time step should be chosen according to the 
duration of the simulation and the precision needed for the spin 
of the inner planet of the system. 

The semi-major axis of the planets shows perfect agreement, 
while the eccentricity, the obliquity, the rotation period, and to 
a lesser extent, the inclination, show a small difference in oscil¬ 
lation frequency. Figure 8 shows the evolution of the obliquity 
of the inner planet of the system corresponding to Case (6”) of 
Table 5 compared with the CR13 code. The mean value of the 



Fig. 7. Evolution of (P p - P p ,o)/F p ,o of the inner planet of the system 
corresponding to Case (6) of Table 5. The colored lines correspond to 
the results of Mercury-T : red for a time step of 0.08 day, green for a 
time step of 0.05 day, light blue for a time step of 0.01, and dark blue 
for a time step of 0.001 day. 


obliquity, as well as the maximum and minimum values, are re¬ 
produced well. The only difference is the oscillation frequency. 

The remaining small difference between Mercury-T and the 
CR13 code, i.e., the rotation period drift and the difference in 
the frequency of the oscillations of quantities, is a result of the 
different integration schemes and numerical effects. 

The system is also very sensitive to the initial conditions. For 
example, for a similar rotation period and obliquity, changing the 
initial direction of the spin of the planet O p leads to a different 
mean value of the oscillations of the obliquity. 

For Case (7), for which all the effects are considered, the 
general behavior is perfectly reproduced. Both planets migrate 
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Time (yr) 


Fig. 8. Evolution of the obliquity of the inner planet of the system cor¬ 
responding to Case (6”) of Table 5 for the last 600 yr of its evolution. 
The black line corresponds to the results of the CR13 code. The red line 
corresponds to the results of Mercury-T. 


outward due to the BD tide and enter a mean motion resonance 
at approximately the same time. By entering the resonance, the 
eccentricities and obliquities evolve similarly. 

We also tested the conservation of energy and total angular 
momentum for these examples. For Case (6), we find that each 
component of the total angular momentum is conserved, and a 
reaches a few 10~ 7 after 10 5 years of evolution. The total energy 
of the system is conserved up to a few 10 -5 . For all the other 
cases, i.e., (6’) to (6’”), the conservation is slightly better, but the 
orders of magnitude are the same. In the case of a non-dissipative 
force, our code conserves total energy and angular momentum. 

Apart from very small differences, Mercury-T and the CR13 
code give the same results. We, therefore, consider this agree¬ 
ment good enough for the study of exoplanets. 

4. The case of Kepler-62 

Just as Mercury-T can be used to study hypothetical systems, 
it can also be used for known exoplanet systems. This code 
has been used in the following articles: Bolmont et al. (2013); 
Quintana et al. (2014); Bolmont et al. (2014b), and Heller et al. 
(2014). Here, we present a study of Kepler-62, a system that 
hosts five planets (Borucki et al. 2013). Two of these planets are 
in the insolation habitable zone (HZ). 

Table 6. Stellar properties 


Mass 

(M 0 ) 

Radius 

(Rq) 

k 2 ,* cr 

P+, 0 

(day) 

0.69 

0.63 

0.03 <x* 

79.7 


Planetary climate depends on many different parameters, in¬ 
cluding orbital distance, eccentricity, obliquity, rotation period, 
and tidal heating (Milankovitch 1941; Spiegel et al. 2009, 2010; 
Dressing et al. 2010). Because all of these parameters are in¬ 
fluenced by tidal interactions, it is of paramount importance to 
consider tides in climate studies (Bolmont et al. 2014a). 

Below, we present a dynamical study of the Kepler-62 sys¬ 
tem, showing the influence of the diverse physical effects, e.g., 
tides, rotation flattening, and general relativity, on the stability 
of the system and on the spin state evolution of the planets. 


To re-create the initial conditions, we used the data from 
Borucki et al. (2013) for semi-major axis, eccentricity, longitude 
of periapsis, and epoch of mid transit. We used the values given 
for the radius of the planets, and we tested the system with dif¬ 
ferent masses for the planets (all of which we assumed to be 
rocky). 

The simulations we show here are, of course, possible evo¬ 
lutions of the system, and we are aware that there are many un¬ 
certainties on many parameters, starting with the masses of the 
planets and their dissipation factors. However, despite these un¬ 
certainties, we aim here to show that some general behaviors can 
be identified in the dynamics of the system. 

4.1. Dynamics and stabilization 

To investigate the effects of tides, rotation-flattening, and general 
relativity on the dynamics of the system, we tested the system for 
five different cases that are listed in Table 8. 

Table 8. Test simulations for stability 


Effects considered 



GR 

Rot. flat. 

Tides 

1 

X 

X 

X 

2 

/ 

X 

X 

3 

/ 

X 

/ 

4 

/ 

/ 

X 

5 

/ 

/ 

/ 


Assuming the planets have a rocky composition (in this hy¬ 
pothesis, planet d has the maximum mass given in Borucki et al. 
2013), we find that the system, hereafter called system J[, is un¬ 
stable in case (1). After 3 Myr, planet c is ejected from the sys¬ 
tem. In this case, planet d is very massive: 14 M®, and its influ¬ 
ence on the less massive planet c destabilizes the system. We find 
that system is also unstable in case (2). However, the desta¬ 
bilization occurs much later, after ~ 20 Myr of evolution. Here, 
the correction for general relativity has the effect of stabilizing 
the system. General relativity also causes apsidal advance, which 
can, therefore lead to situations that are favorable or unfavorable 
to stability, depending on initial conditions. Changing the ini¬ 
tial orbital angles, such as the longitude of periastron, modifies 
the amplitude of the eccentricity oscillations, and this could also 
lead to a more or less stable system. However, for our particular 
choice of initial conditions, it would seem that general relativity 
has a stabilizing effect. 

Using the same masses for planets b, c, e, and f as in system 
Ji, we tested the stability of the system in Cases (1) and (2) for 
different masses of planet d. We found that the system is sys¬ 
tematically stable for masses lower than ~ 8 M® for planet d. 
However for masses higher than ~ 8 M®, most simulations lead 
to destabilization within 30 Myr. For simulations done with a 
mass higher than 8 M® for planet d, destabilization occurs either 
for Case (1) or (2), or both, illustrating the importance of taking 
the correction for general relativity into account. 

We also observed that changing the masses of the planet 
very slightly influences the stability of the system. Changing the 
masses by 5% leads to a stable system in all cases - hereafter 
called system !B. A broad study of the stability of this system is 
beyond the scope of this paper, which illustrates possible future 
studies using the code Mercury-T. The stability should be tested 
for all possible masses, which leads a high number of combi- 
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Table 7. Planetary physical parameters 




Kepler-62b 

Kepler-62c 

Kepler-62d 

Kepler-62e 

Kepler-62f 

Masses (M©) 

3K 

2.60 

0.130 

14 

6.100 

3.500 


S 

2.72 

0.136 

14 

6.324 

3.648 

a (AU) 


0.0553 

0.0929 

0.12 

0.427 

0.718 

ecc 


0.071 

0.187 

0.095 

0.13 

0.094 

inc 


0.8 

0.3 

0.3 

0.02 

0.1 

P p ,0 (hr) 


24 

20 

30 

24 

24 

e P ,o (rad) 


0.1 

0.02 

0.05 

0.03 

0.4 


nations and simulations to have a map of the stability of the 
system (e.g., Laskar 1990; Correia et al. 2005; Couetdic et al. 
2010; Mahajan & Wu 2014). Unstable regions would, therefore, 
correspond to unrealistic configurations. This illustrates the im¬ 
portance of constraining the masses of planetary systems (e.g., 
with HARPS in Dumusque et al. 2014). 

When we add tides, as in Case (3), and assume nominal dis¬ 
sipation factors for the planets, we find that system fft becomes 
stable for the duration of the simulation. In this case, tides have a 
stabilizing effect on the system. Indeed, in our simulations, both 
planetary tides and stellar tides, have a damping effect on the ec¬ 
centricity, therefore reducing the probability of the system to be 
chaotic. System !B remains stable, at least for the duration of the 
simulation of 30 Myr. Semi-major axes and eccentricities do not 
significantly evolve tidally during the simulation, however the 
obliquities and rotation period of the planets do (see Figure 9). 

We tested the evolution of system !B for different planetary 
dissipation factors. The higher the dissipation of planet j, the 
faster its tidal evolution. This effect is first visible on the rota¬ 
tion of the planets (obliquity and rotation period), which evolve 
much more rapidly. It also has a small effect on the eccentricity 
of the planets, which is not visible on the graphs. To quantify 
this, we computed the mean angular momentum deficit (AMD, 
e.g. Laskar 1997) for a set of simulations. 

We did a test, varying the dissipation of each planet from 
a reference simulation corresponding to a dissipation factor of 
0.1 <x© for all planets. We increased the dissipation of each planet 
one after the other from the reference value of 0.1 <x© to 10 <x© 
and 100 cr©. 

The results indicate that, by increasing the dissipation of a 
planet, the AMD decrease. The decrease is more or less pro¬ 
nounced depending on the planet considered. When increasing 
the dissipation of planets c, e, and f to 100 cr©, the AMD is 
~ 0.03% lower than the reference case. However, when increas¬ 
ing the dissipation of planet d, the AMD is 0.055% lower. As 
planet d is very massive in the system, damping its eccentricity 
slightly has consequences for the whole system. When increas¬ 
ing the dissipation of planet b to 100 cr©, the AMD is ~ 0.7% 
lower than the reference case. Increasing the dissipation of the 
closest planet has the biggest effect on the dynamics of the sys¬ 
tem. Increasing the dissipation leads to a slightly less chaotic 
system. 

When we add the effect of rotation flattening, as in Case (4), 
we find that system fft is stable, at least for the duration of the 
simulation of 30 Myr. This effect stabilizes the system by chang¬ 
ing the precession rates. The dynamics of system !B are not sig¬ 
nificantly changed by this effect. 

When we consider all effects , as in case (5), we find that 
System J\ is stable at least for the duration of the simulation 


of 30 Myr. Compared to Case (3), the addition of the rotation¬ 
flattening effect only slightly changes the equilibrium values of 
the obliquities of the planets. 

4.2. Obliquity and rotation period 

Because of the planetary tide, the obliquity decreases and the 
rotation period evolves toward pseudo-synchronization. The 
evolution timescales of these two quantities are shorter than 
the timescales of evolution of semi-major axis and eccentric¬ 
ity. Figure 9 shows that, for Kepler-62, the rotation period of 
the three inner planets of the system evolves toward pseudo¬ 
synchronization in less than 10 Myr, and their obliquities evolve 
toward small equilibrium values (< 1°). 


10.000 


1.000 


13 

cr 




Time (years) 


Fig. 9. Tidal evolution of the Kepler-62 system (S). Top panel: evo¬ 
lution of the obliquities of the five planets. Bottom panel: evolution of 
their rotation periods in solid colored lines. The dashed lines correspond 
to the pseudo-synchronous rotation period and the solid black line cor¬ 
responds to the rotation period of the star. 


These simulations were performed using the values in 
Borucki et al. (2013) as initial conditions, so that our simulation 
shows how the system could evolve in the future. However, we 
can draw some conclusions on the past evolution of the system 
from Figure 9, or from a simple evolution timescale calculation. 
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Given that the age of the system is estimated at 7 Gyr (Borucki 
et al. 2013), we indeed expect that the three inner planets of the 
Kepler-62 system are now rotating slowly (their period is higher 
than 100 hr), and they have quasi null obliquities. 

The ratio between the pseudo-synchronization rate and the 
orbital frequency depends only on the eccentricity of the planet. 
But in a multiplanet system, the eccentricity of a planet is ex¬ 
cited owing to the planet-planet interactions, and oscillates with 
a combination of frequencies that correspond to secular modes 
(e.g., Murray & Dermott 1999). In our simulations, the plan¬ 
ets experience relatively large eccentricity oscillations (see next 
section), causing the planets’ pseudo-synchronization periods to 
also oscillate. In reality, the rotation periods of the planets are 
not exactly equal to the corresponding pseudo-synchronization 
period. Figure 10 shows the evolution of the rotation period of 
Kepler-62b compared to the pseudo-synchronization period and 
the synchronization period. The pseudo-synchronization period 
oscillates too fast for the rotation period to be able to follow. As 
a consequence, the instantaneous rotation period of Kepler-62b 
oscillates out of phase with the pseudo-synchronization period 
and with a lower amplitude. 



9.90x10® 9.92x10® 9.94x10® 9.96x10® 9.98x10® I.OOxlO 7 

Time (years) 

Fig. 10. Short-term (100 000-year) evolution of the rotation period 
of Kepler-62b ( [!B ). Solid line: rotation period. Dashed line: pseudo¬ 
synchronization period. Dashed-dotted line: synchronization period. 


During the 30 Myr of the simulation, the obliquities and ro¬ 
tation periods of the HZ planets Kepler-62e and f did not evolve 
significantly so we performed longer simulations for these two 
outer planets. Assuming an Earth-like dissipation for the two 
planets, we found that Kepler-62e is likely today to have reached 
pseudo-synchronization and have low obliquity. Figure 11 shows 
that after 3 Gyr of evolution, the obliquity has been damped and 
the rotation pseudo-synchronized. For Kepler-62f the timescales 
of evolution are higher and Figure 11 shows that the rotation pe¬ 
riod is still evolving towards pseudo-synchronization after 7 Gyr 
of evolution, and that the obliquity could still be high. 

The dissipation of the planets is not constrained, and chang¬ 
ing the dissipation would only shift the curves right (if the dissi¬ 
pation is lower) or left (if the dissipation is higher). As an Earth¬ 
like dissipation is probably a high dissipation value (due to the 
presence of oceans, e.g. Lambeck 1977), it seems likely that the 
curves should be shifted to the right. 


4.3. Consequence of dynamics on the potential habitability of 
Kepler-62e and Kepler-62f 

Kepler-62e and Kepler-62f are both inside the HZ. However, be¬ 
ing in the HZ does not entail the presence of surface liquid water. 
Surface conditions that are compatible with liquid water depend, 
of course, not only on the properties of the atmosphere (e.g., 
pressure, temperature and chemical composition) but also on or¬ 
bital parameters. These include semi-major axis and eccentricity, 
as well as physical parameters such as the obliquity and the ro¬ 
tation period of the planet (Milankovitch 1941). 

Using the Mercury-T code, we can simulate the dynamical 
evolution of habitable planets within their system, taking into ac¬ 
count the rich dynamics occurring in a multiplanet system. This 
allows us to provide a set of orbital and physical input param¬ 
eters that are consistent with the real dynamics of a planetary 
system, for any kind of climate model. 



0.00 _I_I_I_:_I_I_I_ i_ I _ I _ I _i _I_I_I_:_L 


0 2x10 4 4x10 4 6x10 4 8x10 4 1x10® 

Time (years) 

Fig. 12. Short-term (100 000 year) evolution of the eccentricity of the 
Kepler-62 system planets. 

Figure 12 shows the evolution of the eccentricity of the plan¬ 
ets for 100 000 yr. The eccentricities of Kepler-62e and f oscil¬ 
late respectively between 0.02 and 0.16 and between 0.05 and 
0.19 with a modulated frequency. These important periodical 
changes in eccentricity have an effect on the climate on Kepler- 
62e and f, similar to how the Milankovitch cycles had an impact 
on the paleoclimate of Earth (Berger et al. 1992). 

Furthermore, as seen in the previous section, we can also 
draw some conclusions about the rotation of the planets. We have 
shown that it is likely that Kepler-62e has a pseudo-synchronous 
rotation (or a near pseudo-synchronous rotation, see section 4.2) 
and that its obliquity is very small. This kind of planet, with a 
slow rotation (almost 3 000 hour or 125 day), could have large 
Hadley cells that bring hot air to the poles and there would be no 
longitudinal circulation (Merlis & Schneider 2010; Leconte et al. 
2013). However, as Kepler-62e is close to the inner boundary of 
the HZ, the surface temperatures might not reach low enough 
values to create cold traps. Consequently, a study of this planet 
using a global circulation model would be needed to test the po¬ 
tential of this planet to host surface liquid water. 

Figure 11 shows that Kepler-62f could have a high obliquity 
and a fast rotation period (as fast as the Earth’s 24-hour rotation). 
Of course, we do not know the initial conditions on the spin of 
the planets. However, formation scenarios show that planets are 
likely to have an initial fast rotation rate because of collisions 
and an isotropic distribution of obliquities (Kokubo & Ida 2007). 
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Kepler-62e Kepler-62f 



Bottom panel: evolution of their rotation periods in colored solid lines. The dashed lines correspond to the pseudo-synchronous rotation period, 
and the solid black line corresponds to the corotation radius. 


For fast rotation rates, the obliquity is excited and can reach high 
values even if it started at a low value (see, for example, the solid 
red line in Figure 11). As a result, there is therefore a high proba¬ 
bility that the obliquity of Kepler-62f is actually non-negligible. 
Furthermore, its rotation period could still be quite fast: between 
20 and 40 hr at the assumed age of the system (see Figure 11). 
It is therefore likely that Kepler-62f would have a very different 
type of climate from its neighbor. Indeed, a non-negligible obliq¬ 
uity would lead to seasonal effects, and a fast rotation would lead 
to a different wind pattern with not only latitudinal winds, but 
also longitudinal winds. 

5. Conclusions 

In this study, we have presented a code that computes orbital 
evolution for tidally evolving multiplanet systems. The theory on 
which this code is based is the constant time lag model, which is 
an equilibrium tide model. This code allows the user to compute 
the evolution of the orbital distance, eccentricity and inclination 
of planets, as well as their rotation state (obliquity and rotation 
period). It also computes the rotation period of the host star con¬ 
sistently (taking into account the spin-up due to radius-shrinking 
and the effects of tides). 

The evolution tracks of the radius of various host bodies were 
implemented: BDs of mass between 0.01 and 0.08 M©, M dwarfs 
of 0.1 M©, Sun-like stars and Jupiter. This allows the user to 
study the influence of a changing radius of the host body on the 
tidal evolution of planets. 

In this work, we have endeavored to validate our code. To 
this end, we compared the outputs of a code that solves the tidal 
secular equations of single-planet systems (see Bolmont et al. 


2011, 2012), with the outputs of our new code. We also tested 
the rotational-flattening effect, by comparing Mercury-T with the 
CR13 code, which was developed independently. We found that 
Mercury-T reproduces the secular evolution of the planets well. 
We also made sure that the total angular momentum was con¬ 
served in all of our examples. 

Potential users of this code should bear in mind that, when 
a planet is alone in the system, the code can produce a spuri¬ 
ous remnant eccentricity. For each simulation, we also advise the 
user to verify the conservation of total angular momentum and 
the robustness of the spin integration by doing a simulation with¬ 
out tides and with the effect of the rotational-induced flattening 
to see whether there is a drift of the mean value of the obliquity. 
If there is a drift, then the time step has to be decreased. 

Some ongoing improvements in this code would consist of 
improving the models of the planets. Indeed, the use of the 
constant time lag model for terrestrial planet has been criti¬ 
cized (e.g., Makarov & Efroimsky 2013; Efroimsky & Makarov 
2013; Makarov & Berghea 2014; Correia et al. 2014), and it 
is probable that the planets are not evolving towards pseudo¬ 
synchronization, but are trapped in spin-orbit resonances. Be¬ 
sides, to correctly determine the spin of planets, one needs to 
take thermal tides into account (e.g. Cunha et al. 2014; Leconte 
et al. 2015). With a global circulation model, Leconte et al. 
(2015) showed that this phenomenon can drive planets out of 
synchronization even if they have a thin atmosphere. We also in¬ 
tend to implement a better description of the dissipation within 
the star (e.g., using models found in Auclair-Desrotour et al. 
2014). A wind prescription will also be added soon (as in Bol¬ 
mont et al. 2012). In the future, we also intend to investigate the 
multibulge effect, i.e., the influence of the bulge raised on the 
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star by planet j on the dynamical evolution of planet i Aj (as was 
achieved in Touma & Wisdom 1994, for the Earth-Moon-Sun 
system). 

Mercury-T is a very powerful tool for simulating the evolu¬ 
tion of any kind of planetary system. It can be used to simulate 
known exoplanetary systems to try to: identify trends, as we have 
done in this work for the Kepler-62 system; investigate the stabil¬ 
ity of the system, taking all the important physical phenomena 
into account; and to investigate the influence of tidal dissipa¬ 
tion factors on the evolution of the system, to maybe constrain 
the parameters space. For example, Bolmont et al. (2013) used a 
previous version of the code to evaluate the possible eccentric¬ 
ity of the transiting inner planet of the 55 Cancri system. Us¬ 
ing this code, they investigated if tidal heating could contribute 
significantly to 55 Cancri e’s thermal emission. Our code also 
allows the user to have an idea of the spin state of planets (as 
in this work or in Bolmont et al. 2014b, which focuses on the 
Kepler-186 system and also constitutes a fine example of the use 
of Mercury-T). 

Mercury-T is particularly interesting to use for simulating 
the orbital dynamical evolution of habitable planets because it 
allows for reasonable and consistent input in climate models to 
investigate the potential of these planets to host surface liquid 
water and also to investigate the influence of eccentricity oscil¬ 
lations on such a climate. 
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